Local minimal energy landscapes in river networks 
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The existence and stability of the universality class as- 
sociated to local minimal energy landscapes is investigated. 
Using extensive numerical simulations, we first study the de- 
pendence on a parameter 7 of a partial differential equation 
which was proposed to describe the evolution of a rugged land- 
scape toward a local minimum of the dissipated energy. We 
then compare the results with those obtained by an evolution 
scheme based on a variational principle (the optimal channel 
networks). It is found that both models yield qualitatively 
similar river patterns and similar dependence on 7. The ag- 
gregation mechanism is however strongly dependent on the 
value of 7. A careful analysis suggests that scaling behaviors 
may weakly depend both on 7 and on initial condition, but 
in all cases it is within observational data predictions. Con- 
sequences of our results are finally discussed and the most 
plausible scenario is presented. 

92.40.Fb,64.60.Ht,05.60.+w 



I. INTRODUCTION 

Understanding the development of a landscape 
in terms of fundamental mechanical principles is a 
formidable task In spite of this high complexity, re- 
cent theoretical studies have resulted in considerable pro- 
gresses by considering the issue from a viewpoint anal- 
ogous to the one taken in conventional critical phenom- 
ena where simple models are exploited to identify 
universality classes. The main idea of this approach is 
indeed to focus on few fundamental ingredients which, in 
the spirit of critical phenomena, are expected to provide 
a reasonable description of the large-scales properties of 
the network. 

The remarkable properties of river networks have been 
known for some time and they can be condensed in 
few phenomenological scaling laws which have been con- 
firmed to hold in the observational data ||. While 
these laws do not explain the underlying physical mech- 
anisms, they nevertheless provide guidelines for their 
search. Hence any newly proposed model for river net- 
works ought to be tested against these laws. In the lan- 
guage of critical phenomena, those scaling laws can be 
used to derive critical exponents and thus discriminate 
among different universality classes. Among such laws a 
fundamental role in the physics of the erosion of a land- 
scape due to the flow of water over it, is played by the 
slope-area law (see Ref. jl] for a review) . 



A few Langevin-like equations have been recently 
proposed for describing the evolution of the landscape 
under the effect of erosional processes. In Ref. D, 
reparametrization invariance arguments Jl0| were used 
to derive a dynamical equation which yields the slope- 
area law as a stationary state. The same equation was 
obtained by Somafai and Sander using Landau ar- 
guments. Other proposals have also been advanced ||, 

The equation proposed in Ref. H, was studied both 
analytically in 1+1 dimensions (one spatial and one tem- 
poral) and numerically using a self-consistent (SC) solu- 
tion in 2+1 dimensions, and it was shown to predict a 
fairly reasonable stationary state quite different from the 
starting network acting as initial condition. This type of 
analysis hinges on the observation that the network rep- 
resenting the river has an evolution mirroring that of the 
evolution profile. Other methods have also been devised 
which directly address the topological properties of the 
network itself, the best known being the optimal channel 
networks (OCN) This is a lattice model where a 

functional describing the dissipated energy is minimized 
in order to find the optimal configuration and it is based 
on the idea that, presumably, the erosional process tak- 
ing place on a landscape is driven by a strive for opti- 
mality. A simulating annealing procedure Jig] has been 
implemented to this aim both at zero and finite tem- 
perature |16|. While the latter is aimed to find the ab- 
solute minimum, the former is expected to display local 
minima only. Using exact bounds and finite-size scaling, 
Maritan et al ( |T^,|l8|) showed that the absolute mini- 
mum belongs to the mean-field universality class, that, in 
turn, means that the corresponding network has a highly 
symmetric pattern with small rivers draining into bigger 
rivers in a predictable way (this network is akin to the 
Peano basin |19|). However this minimum is not eas- 
ily reachable in the space of all configurations, and one is 
then led to suspect that real rivers are better described by 
configurations related to local minima. We further note 
that stationary states of the aforementioned dynamical 
equations, are also expected to be associated with local 
minima for the same reason. 

In view of the above discussion, it is apparent that a 
more complete description of the stability and the scaling 
behavior associated with these local minima would be de- 
sirable. The present work is an attempt in this direction. 
Our aims are threefold. Our first goal is the full char- 
acterization of the possible universality class associated 
with local minima. Using the SC numerical procedure 
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envisag ed in Ref. @, we extend previous numerical re- 
sults both in size and statistics and we use a finite size 
procedure for a more accurate estimate of the exponents. 
We then compare these exponents with those calculated 
using a zero-temperature simulation in the OCN frame- 
work with similar sizes and statistics. Our second goal is 
the generalization of the dynamical equation to include 
a tunable parameter which turns out to be related to the 
< 7 < 1 parameter considered in OCNs and which is 
responsible of the aggregation mechanism. In Ref. fl§| l 
it was shown that the absolute minimum is insensitive 
to a variation of 7 in the interval 1/2 < 7 < 1 [ p0[ . 
We find that although the final patterns display marked 
differences as a function of 7, critical exponents show a 
much smaller discrepancy, which our results indicate to 
be marginal (at the edge of error bars), both in the SC 
and OCN case. Finally we test the stability of our results 
against a variation in the initial conditions. We find that 
even in this case a marginal difference appears in the 
effective critical exponents. 

The paper is organized as follows. In Secjnl, we re- 
view the dynamical equation whereas in Sec.pl scaling 
laws and critical exponents are briefly recalled. Sec. IV 
contains the definition of OCN and our results are pre- 
sented in Sec.[y|. Finally Sec. VI contains some concluding 
remarks and future perspectives. 



II. THE DYNAMICAL EQUATION 

A rather general dynamical equation consistent with 
general principles, and capturing the physics of erosional 
processes occurring in river basin, is |7j]: 



\Vh(x,t)\ -Q^OM) 



(2) 



dh(x,t) 
dt 



= u + vV 2 h(x,t) - Q a (x,t) Vh(x,t) +n(x,t) 



Here h(x,t) and Q(x,t) are the height of the landscape 
and the flow rate at position x and time t, r)(x, t) is a noise 
term and a is a parameter which will be related to 7 (see 
below). The constant term u on the right-hand-side of 
Eq. (||) mimics the so-called geological uplift which is 
known to originate by tectonic forces. The second term 
represents a local diffusion term of strength v, condensing 
both smoothing (rounding of hilltops) and sedimentation 
(filling of valleys) processes. The third term is non-local 
and corresponds to an erosion driven by the waterflow on 
a surface with Q(x,t) representing the flow through site 
x. The exponent a is a parameter of the model which was 
assumed to be equal to 1 in Ref. || and Ref. j?J . The last 
noise term is added to account for small-scale stochastic 
processes such as rainfall fluctuations. In the attempt 
of extracting the basic features associated with the third 
term, one can ignore both the diffusion and noise terms 
(the accuracy of this approximation has been discussed 
111 Ref. j§ ). In this case, a stationary state is obtained 
when the uplift term balances the flow dependent ero- 
sional term, that is when 



where, for later convenience, we define 7 = 1 — a/2. Most 
of previous work was performed with a = 1 (7 = 1/2) 
(there are few exceptions but with aims and methodolo- 
gies different from ours, see e.g. Ref. 0). In this case Eq. 
(||) is known as slope-area law Jjj , which is a robust em- 
pirical law verified by real field observational data. We 
shall see later on how a is related to the dissipated en- 
ergy functional. Here we note that if we assume (as we 
shall do hereafter) uniform rainfall (and no ground wa- 
ter), then Q{x, t) ~ a(x, t) where a(x, t) is the area of the 
basin draining into point x at time t. The basic result 
of this dynamical equation is that the non-linear term is 
able to account for the correct relationship between local 
slope and water flowrate. 



III. CRITICAL EXPONENTS AND SCALING 
LAWS 

River networks are a remarkable example of systems 
governed by scaling laws. Although the concept of power 
and scaling laws has been known to the hydrologists for 
half a century, it was not until recently that this con- 
cept was put into a well defined framework 0, |Q in 
analogy with conventional non-equilibrium critical phe- 
nomena where power laws are associated with critical 
exponents. For the sake of completeness, we shall briefly 
review here a few of the central laws appearing in river 
networks which we regards as the most fundamental. 

Hack observed that the total area a draining into a 
given point and the upstream length I going from that 
point to the source though the path of maximum water 



^jlow, were not independent but related by j21| 

I ~ a h 



(3) 



where h is often referred as Hack's exponent and it ranges 
in the interval [0.5 — 0.6] in real rivers. The distribu- 
tions of drainage areas and upstream lengths also follow 
a power law. Within the context of finite-size scaling, 
these may be written as 



for areas and 



P ( ai L)=a- T f(a/L*) 



n(l,L) = r^g(l/L d ') 



(4) 



(5) 



for lengths. In Eq. (|J) p(a, L) is the distribution of 
drainage areas on a basin of size L and f(x) is a finite size 
function. It defines two critical exponents r and <j> which 
are not independent Q. Similarly n(l, L) is the distribu- 
tion of the upstream lengths on a basin of size L, g(x) is 
a finite size function and there exists a relation between 
ip and di. Note that (j> and all define the "typical" area 
(otyp ~ L^) and length (l typ ~ L dl ). Finally we remark 
that one usually distinguishes between self-affine basins 
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(di = 1, cj) < 2) and self-similar basins (di > 1, = 2) 
|B2j and in both cases only one exponent out of four is 
independent. They are related via scaling relations as 



1 
h 



(G) 



IV. 



ZERO-TEMPERATURE OCN: LOCAL 
MINIMA 



Within the context of the OCNs, the "dissipated en- 
ergy" of a river network at a given time, can be written, 
in continuum notations, as (Til: 



(7) 



E(t) = / dx \Vh(x,t)\Q(x,t) 



The local gradient |V/i(x,i)| and the water flow rate 
Q(x,t) are expected to satisfy Eq. (0). Hence one gets: 



E(t) 



dx Q 1 {x, t) 



(8) 



We note that although the energy expression (||) has been 
derived in the context of OCN in principle it could 
be defined in the evolution equation ([!]) as well, thus 
providing a way of monitoring the rate of approach to 
steady state. 

The basic assumption of the OCN is that there is a 
tendency of any real river basin, to assume a configura- 
tion which minimizes Eq. (||). A natural question arising 
is the characterization of the local and absolute minima 
of E. ft was shown JT^Jiql , that the absolute minimum 
of E is insensitive to a variation of the value of 7 in the 
interval 1/2 < 7 < 1 [§(]] and then h = 1/2, r = 3/2 and 
di = 1 . The analogous issue in the case of local minima 
is much less clear. In fact, when 7 = 0.5, numerical work 
fl8f suggests that there exists a set of local minima which 
presumably corresponds to a new universality class that 
is the relevant one for real rivers. Indeed only exceptional 
events are able to radically modify river courses and so 
local minima of dissipated energy trap the system with 
high probability and therefore dominate the statistics. 

It is then a vital issue to discriminate whether or not 
local minima configurations are indeed related to a well 
defined and robust universality class independent on 7. 



V. RESULTS 

In this section, we describe numerical procedures and 
results in detail for each case. All calculations are car- 
ried out on a L x L square lattice with periodic bound- 
ary conditions in one direction which we identify as the 
transverse direction. Multiple outlets are allowed in the 



outflowing longitudinal direction (West side in the fig- 
ures) whereas an infinite wall is set up on the opposite 
side (East side). All averages considered in the statis- 
tics are carried out only over the river with the largest 
flow. It is worth stressing that the choice of considering 
only the maximum river in a multiple outlets environ- 
ment corresponds to consider the statistical behavior of 
rivers that are in competition to drain a given region. 
On the other hand, a single river within a given region is 
more appropriate if geological constraints are known to 
exist. 

Typically 8 nearest-neighbors (nn) - 4 associated with 
the square lattice and 4 associated with the two diago- 
nal directions - are allowed. A somewhat more restricted 
choice considers only the 4 natural nn associated with the 
square lattice structure. In view of universality one would 
expect that the details of the lattice structure should not 
matter. This second choice has the considerable advan- 
tage of being less time consuming for numerical purposes. 
For this reason, although all following figures display pat- 
terns obtained with 8 nn, the results reported in this 
work are obtained with statistics based on 4 nn. In one 
example, we have explicitly checked that the outcomes 
using the two choices are consistent within the statistical 
errors. 

Finally, in all our networks, the drainage area a(x, t) is 
computed, at each time step, in a standard way according 
to i: 



a(x,t) = ^2a(y,t) + 1 



(9) 



y(x) 



where y(x) denotes all sites y which drain into x and the 
last term on the right hand side represents a unit rainfall 
input on each site at all times. 



A. Initial condition 

There are many possible initial conditions that can be 
used in numerical analysis of river networks. A popu- 
lar one among hydrologists is a deterministic comb-like 
structure fl. However this initial condition suffers of 
various drawbacks [ p3| and its use in the multiple outlets 
case, such the one treated here, appears to be somewhat 
inconvenient. Hence we consider here two other phys- 
ically reasonable choices, namely a spanning tree (ST) 
and a Scheidegger network (SN). Although STs are well 
known in statistical physics mainly due to their relation 
with the q — > limit of the Potts model pij , only recently 
their topological properties have been studied in details 
|p5| . A suitable variation of spanning trees has even been 
proposed as a topological model for river networks p2[ , 
fl26| . We have generated STs with multiple outlets by us- 
ing an adapted Broder's algorithm (see Ref. |2^| and ref- 
erences therein). We have considered sizes ranging from 
L = 32 to L = 512 and statistics based on a number of 
configurations ranging from 500 to 100 respectively. We 
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have then performed a finite size analysis described in 
the next subsection to extract the most reliable values 
of the exponents. Our best result for the exponents are 
t = 1.378 ± 0.002, V = 1-596 ± 0.003, h = 0.633 ± 0.003, 
di = 1.25 ± 0.01, (f> = 2.00 ± 0.01 in excellent agreement 
with the exact results r = 1.375, ip = 1.6, h = 0.625, 
di = 1.25 p5fl . This also provides a good test on the 
quality of our data analysis. Note that the exact value 
of 0, albeit not known, can be derived from the knowl- 
edge of h using Eq. (||). This predicts = 2.0 in perfect 
agreement with our computed value. We also note that 
the obtained results are nearly identical to the one pre- 
dicted and numerically generated for a single outlet [^6| . 

A ST is a self-similar network in that it is undirected 
and isotropic. A quite different choice (directed and 
anisotropic and hence self-affinc) is a SN. This was again 
proposed as a topological model for river networks (see 
e.g. [jjj) and the r exponent has been exactly determined 
via a mapping to a one-dimensional model for mass ag- 
gregation pTij ]. A similar mapping to a diffusion-reaction 
model also provides the solution in general dimension- 
ality ]2^]. Our best results for the critical exponents 
with the same statistics as above are r = 1.337 ± 0.003, 
ip = 1.49 ± 0.02, h = 0.69 ± 0.02, di = 0.96 ± 0.01, 
4> = 1.53 ± 0.02, which are in excellent agreement with 
the expected values (r being exact, the others determined 
via Eq. (§)), that is r = 4/3, ip = 1.5, h = 0.75, d l = 1, 
4> = 1.5. 

In Figs|l] and |^, we depict typical patterns for a ST 
and a SN respectively. The presence of multiple outlets 
magnifies the main difference between ST and SN. A ST 
is constructed in such a way that total freedom is given 
to the meandering of streams, the only constraint being 
that they all terminate on the same line (West side in 
the Figures). Typically this allows the formation of a 
main big river of size considerably larger with respect to 
the others. For SNs, the East- West preferred direction 
prevents the forming of such big river and many smaller 
rivers are usually present as shown in Fig.|[ 



B. The Dynamical Equation: A Self-Consistent 
solution 

As we discussed in Sec. [n|, we seek the stationary 
states of the following simplified equation: 



dh(x, t) 
dt 



= u-Q a (x,t) Vh(x,t) +r)(x,t) 



(10) 



where a = 2(1 — 7). The stationary averaged states of 
Eq. (|l0|) are expected to conform to Eq. (g). Despite the 
apparent simplicity of the equation, an explicit numerical 
solution proves to be rather slow The reason for 
this can be traced back to the particular form of the 
erosion term. According to Eq. ( [To| ) only sites with a non- 
negligible combination Q a {x, t) \Vh(x, t)] 2 will affect the 
change of the pattern. As it turns out, this yields a long 



time transient during which the elevation profile evolves 
very slowly. 

Since we are mainly interested in stationary states, 
there is a way out from this situation ||. The main 
idea is to start with an arbitrary network (e.g. a ST or a 
SN) and recursively construct the heights starting from 
the outlets with the aid of Eq. (|^). From the derived 
landscape, a new network (in general different from the 
original one) can then be obtained by assuming that at 
each site the outward direction is along the steepest de- 
scent path and using Eq.(||). The noise term in Eq. ( |l0| ) 
is mimicked by the unity term in Eq. (||). The proce- 
dure can then be iterated until self-consistency is finally 
achieved. The final configuration is, by definition, a sta- 
tionary state of Eq. ([To|). 

The convergence of the above procedure as a function 
of the number of iterations is reported in Fig. [|for various 
values of 7. Two remarks are in order. First we note 
that the dissipated energy as obtained from the above 
SC procedure does not have any physical meaning in the 
transient state. In other words the definition of "time" 
for the independent variable appearing in Fig. |3| should 
be considered only as a short-hand notation for "number 
of iterations". Second it is apparent a the value 7 = 0.5 
seems to be the one with the slowest convergence ratio. 
This is probably due to the particular role played by the 
value 7 = 0.5 as it will become clear shortly. 

Figs.HH depict typical patterns obtained on changing 
the parameter 7 in the interval [0,1]. The effect of the 
value of 7 on the aggregation pattern is evident. As 7 
increases from to 1, single big rivers draining the entire 
basin in a snake-like form are less and less favored. The 
pattern for 7 = has a strong memory of the original 
initial tree. The overall effect of the SC procedure is to 
disfavor long meandering of the streams thus providing a 
sclf-affine character to the final tree. The main river be- 
comes rather straight for 7 = 0.25 and the whole pattern 
appears to be more symmetric with respect to the case 
of 7 = 0. A noteworthy feature is that this tendency is 
inverted as 7 — > 0.5 and returns back for 7 = 0.75. As 
7 — ► 1 rivers become very directed as one would expect 
on the ground that this limit corresponds to the Schei- 
degger network [fl8|| . The fact that 7 = 0.5 most closely 
resembles real rivers is a reflection of the natural selec- 
tion by the erosional processes of this value of 7 in terms 
of Eq. (j^) as it has already been well documented in the 
literature (see e.g. |Q). 

In our numerical estimates of the exponents, for sim- 
plicity, we shall restrict our attention to the half region 
[0, 0.5] . We also note that this is the region inaccessible 
to the analytical scheme of Ref . |17| . 

For a more accurate evaluation of the critical expo- 
nents r and ip it proves convenient to introduce the inte- 
grated probabilities: 

P(a, L) = jf da'p(a', L) = ^F^) (11) 

and 
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m,L) = J o dl'-K{l',L) = l 1 -^G{^-) 



(12) 



where F{x) and G{x) are related to f(x) and g(x) de- 
fined in Eqs. (|) and (§), in an obvious way. An efficient 
way of computing the exponents is through the so-called 
"effective" (sometime also referred to as "running") ex- 
ponents. In the present case they are defined as: 



t(o) = 1 - 



d log P(a,L) 
d a 



and 



d \ogU{l,L) 
dl ' 



(13) 



(14) 



One then obtains an effective exponent for each value of 
the independent variable (a or /). 

Fig. [9] shows one typical result on a 256 x 256 lattice. 
We can divide the obtained values roughly in four regions. 
The first region (1 < a < 10) corresponds to the region 
of no scaling. Small rivers belonging to the second region 
(10 < a < 100) have an exponent close to the absolute 
minimum value r = 1.5. This is consistent with the pic- 
ture of typical rivers (see Fig. |^) where small rivers dis- 
play a marked straightness similar to the one of the abso- 
lute minimum |i"8| and it means that they quickly assume 
configurations consistent with the absolute minimum (in- 
dependent of the initial conditions). Larger areas are as- 
sociated with larger rivers which have longer memory of 
the initial condition (ST in the present case). Hence the 
corresponding exponent is sensibly smaller (closer to the 
ST value 1.38). The last region corresponds to the finite 
size cut-off and must be discarded. 

After discarding the first and forth regions the ob- 
tained values can then be grouped into local bins and 
a local average exponent can be associated with each of 
them. Statistical fluctuations within each box then yield 
an estimate of error bars. This provides our best esti- 
mate of the exponent for each value of L and a simple 
1 / L extrapolation is then carried out to extract the final 
values. This is depicted in Fig. [l^ and the corresponding 
best estimates of this method are reported in Table Q. 
An analogous procedure leads to the best estimates for 
tp as reported again in Table |. Sizes and statistics are 
identical to those considered for the initial conditions and 
hence these simulations are rather time consuming. 

One can notice a weak dependence on 7 for both r and 
ij). One way of computing <j) and di is through the collapse 
plots of the probabilities P(a,L) and H(l,L). However 
we have found that a satisfactory collapse can achieve 
only within a limited range of the appropriate variable 
(a/L* &ndl/L d '). Hence we have opted for an alternative 
scheme hinging on the calculation of the following ratios: 



(- q+l )a 



1,2, 



(15) 



where averages are over the probability densities p(a, L) 
(of the maximum river) and the L dependence is straight- 
forward Bj. A similar relationship holds for the lengths: 



MHL) = 



L° 



1,2,. 



(16) 



The results for q — 1,2,3 are reported in Tabic |Hj for 
the exponent <p and in Table IV for the exponent di . Sur- 
prisingly the values for <fi are consistently larger that the 
space-filling value <j> = 2 which is expected to be the up- 
per bound for this exponent. This is probably due to a 
deficiency of this procedure during a cross-over from a 
self-similar regime (the initial ST) to a self-affine pattern 
(the final configuration). Indeed we shall see later on 
that this feature is not present when one starts with a 
self-affine network (e.g. a SN) from the outset. A seem- 
ingly large value of di is appearing probably due to the 
same reason. Overall, one can notice a rather weak (if 
any) dependence on 7 and almost no dependence on q. 

Finally Hack's exponent h was computed from its def- 
inition Eq. (|^) using an effective exponent method to be 
described below. 

Let us assume a power-law dependence for a generic 
function as given by 



/(*) 



> 



(17) 



On integrating f(x) (between lower and upper limits, say 
xq and x), we find 



xf(x) - x f(x ) 



(18) 



The effective exponent obtained using Eq.(18) in Eq.(|3j) 
with f{x) = a, x = I, and 9 = 1/h, can then be ana- 
lyzed with the same procedure (local average plus 1/L 
extrapolation) outlined before. We note that for a we 
have used an averaged value over all areas correspond- 
ing to the same length. The final results are reported 
in Table Q and one can see that there is an overall good 
agreement with scaling laws. 



C. OCN 

The minimization of the energy functional Eq. (^) goes 
through an algorithm akin to the one exploited by Sun 
et. al ]16[ ]. It is based on the following steps: 

1) An initial configuration (ST or SN) is generated 

and its dissipated energy is computed according to 
Eq. (||). By definition, this initial network has no 
loops. 

2) A link is randomly selected and its local outflow is 

also randomly chosen. 
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3) This new configuration is tested for loop creation. If a 

loop has been created, the configuration is rejected 
and step 2) is repeated. 

4) The energy of the new candidate configuration is com- 

puted. If it is smaller than the previous one, the 
new configuration is accepted, otherwise it is re- 
jected. 

5) Steps 2)-4) are repeated until the energy does not 

change within a given tolerance. 

The final configuration is regarded as a local minimum. 
This scheme is patterned after a standard Metropolis al- 
gorithm at zero temperature, and averages are over many 
different configurations (ranging from 500 at L = 32 to 
100 at L = 256). Fig. [y] depicts the dissipated energy 
per unit of length as a function of the convergence "time" 
(i.e. the number of total iteration of the algorithm) . The 
similarity with patterns obtained from the SC procedure 
is evident. Here, however, the typical convergence time 
is much longer than before. 

Critical exponents are computed with the same pre- 
scription given in the previous subsection. In Fig. [l^ 
we show the resulting 1/L extrapolation. The final best 
estimates are reported in Table |. 

Once more, a weak dependence on 7 (r increases as 7 
increases) can be noticed. All other exponents are consis- 
tent with scaling relations Eq. @. It is worth stressing 
that exponents are nearly consistent (at the edge of sta- 
tistical errors) with those previously obtained from the 
dynamical equation. 

Reg ardi ng the exponents 4> and di , they can be found in 
Tables III and IV respectively. Here too the same feature, 
discussed in V B in connection with the exponent values 
of 4> and di, applies. 



D. Independence of the initial condition 

Our final task is a test of the sensitivity of critical ex- 
ponents to the initial conditions. To this aim we have 
changed the initial condition from a ST to a SN for both 
the dynamical equation and the OCN. The main differ- 
ence between the two initial conditions is that while the 
latter is a directed network, the former is not. Table 
reports the comparison for the case 7 = 0.5, and Fig. Il3 
depicts the network resulting from the SC scheme, for a 
typical final state. Despite the obvious memory of the 
initial network (see Fig. |J), the typical distance among 
big rivers is larger than the initial one. This in turn yields 
a higher value for H and hence non-trivial exponents. 

It is remarkable that although these patterns appear 
to be considerably different from those obtained when 
starting with a ST, the two sets of critical exponents are 
nevertheless consistent within the statistical errors. 

As a final comment, we find in this case values of <p 
and di in very good agreement with scaling predictions. 
The results for different ratios q in the case 7 = 0.5 for 



both the SC and the OCN, along with the comparison 
with the corresponding values stemming from STs are 
reported in Tables [v] and [v| respectively. 



VI. CONCLUSIONS 

In this paper we have addressed the issue of the exis- 
tence and robustness of the universality class associated 
with landscapes corresponding to local minimal energies. 

To this aim we have first extended previous studies 
for both the SC solution of a Langevin equation and the 
OCNs variational methods. Higher sizes and statistics 
have been exploited in both cases and, (to our knowledge) 
for the first time, the most physical procedure of basing 
the statistics only on the river with largest flow, has been 
implemented. 

Second, we have monitored the dependence of critical 
exponents on a parameter 7 associated with the slope- 
area law in one case (SC) and with the dissipated energy 
in the other (OCN). Our results give compatible critical 
exponents between SC and OCN within the error bars, 
but a weak and similar dependence on 7 appears in both 
models. 

Finally we have tested the stability of the obtained 
results for both the SC and the OCN with respect to 
changes in the initial conditions. Although the obtained 
final patterns display a dependence on initial conditions, 
critical exponents appear to be insensitive to this depen- 
dence. 

As a by-product of our investigation, we have found 
that the SC solution of the dynamical equation is a very 
powerful method to investigate river networks as it is ca- 
pable of providing useful informations on the stationary 
state in a simple and physical way. Another interesting 
point, from a numerical point of view, is that this pro- 
cedure typically achieves convergence much faster with 
respect to OCN scheme. 

In view of our results we can now summarize the argu- 
ments favoring and disfavoring the appearance of a new 
universality class associated with local minima. As we 
mentioned in our discussion of Fig. ^ the typical evolu- 
tion of a network appears to depend on the considered 
length scale. Small rivers very quickly settle to their final 
state whereas much longer time is required to large rivers 
to forget their initial conditions. This is also reflected by 
the difficulty in collapsing the distribution probabilities 
of both a and I into a single plot for a reasonably ex- 
tended range of the corresponding variable. It is then 
possible that although only two universality classes (one 
associated with the initial condition and the other to the 
absolute minimum) are present, an intermediate univer- 
sality class sets on due to both the averaging over dif- 
ferent regimes and the difficulty of reaching the absolute 
minimum. 

On the other hand, in support of the existence of a new 
universality class, we cite the fact that critical exponents 



G 



are different from those of both ST and SN - on start- 
ing with these initial conditions, critical exponents in the 
final configurations clearly deviate from their initial val- 
ues. Furthermore all exponents are found to be robust 
and obey to scaling relations summarized in Sec.|T|. 

Overall we believe that the evidence suggested by our 
results favors this second possibility which has also been 
hinted into a different context (£9|. In this respect the 
weak 7 dependence of critical exponents remains unex- 
plained. 

It would be interesting to generalize the present cal- 
culation in two aspects. For the dynamical equation, it 
would be instructive to tackle the problem of the ex- 
plicit solution of Eq. (|l]) . This route has the advantage 
that Eq. (^) (the key relation between erosional process 
and network topology) can be then derived rather than 
assumed as in the self-consistent procedure. Similarly, 
a parallel calculation could also be implemented in the 
OCN framework, upon starting with the more general 
expression given in Eq. (^). 
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FIG. 1. A spanning tree on a 64 x 64 square lattice with 8 
nearest-neighbors, 4 in the Nord-South and East- West direc- 
tion and 4 along the diagonal directions. The outlet is on the 
West side (outflowing), on the East side there is an infinite 
wall, whereas there are periodic boundary conditions in the 
North-South direction. The thickness of the line at each point 
is proportional to the flow through that point. The seeming 
loops are just an artifact of the drawing. 




FIG. 2. A Scheidegger network with L = 64. Boundary 
conditions are the same as in Fig^j. The directness of the 
network provides a privileged East- West direction. 



7 



FIG. 3. Dissipated energy per unit of length, as a function 
of the number of iterations (time t) of the self-consistent solu- 
tion of the dynamical equation, on starting with the spanning 
tree of Fig IlL The parameter 7 is the value appearing in Eq. 
(2). The system size is L = 512 and each point of the curve 
is an average over all configurations which have gone at least 
that far in the number of iterations. All quantities reported 
in this and following figures are dimensionless. 



FIG. 5. Same as in Fig. W, with 7 




FIG. 4. A typical network obtained as a fina l output of the 
self-consistent procedure described in Sec. |VB. 




with 7 



. Here L = 64, 

7 = 0, and the initial condition is the spanning tree of Fig. m 
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FIG. 8. Same as in Fig. W, with 7 
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FIG. 10. Finite size 1/L extrapolation for the value 
of t obtained with the dynamical equation for the cases 
7 = 0, 0.25, 0.5. In all cases the initial condition is the span- 
ning tree of Fig. |l|. 
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FIG. 11. Dissipated energy per unit of length, as a function 
of the number of iterations in the OCN procedure. Here the 
size is L = 256 and again each point of the curves are averaged 
over all configurations which have reached that time t. The 
initial condition is the spanning tree of Fig. [|. 
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FIG. 13. A typical 64 x 64 network obtained as a final 
output of the self-consistent procedure described in Sec. VB 
upon starting with a the Scheidegger network of Fig. 0. 



TABLE I. Critical exponents r, ip and h as a function of 7 
for both the Self-Consistent (SC) procedure and the optimal 
channel networks (OCN) at zero temperature. The values in 
the parenthesis are those obtained from the computed value 
of t and using the scaling relations. For 7 = 0, the OCN 
exponents reported are the exact ones corresponding to the 
initial conditions (a spanning tree in the present case) since 
the energy is a constant. 



7 



TQCN 



TSC 



IpOCN 



IpSC 



0.00 1.38 1.40 ±0.01 1.6 1.69 ± 0.03 (1.67) 

0.25 1.42 ±0.01 1.42 ±0.01 1.71 ± 0.07 (1.72) 1.68 ± 0.05 (1.72) 0.64 ± ( 
0.50 1.44 ±0.01 1.46 ±0.01 1.8 ±0.1 (1.79) 1.82 ±0.05 (1.85) 0.61 ± ( 



TABLE II. Critical exponents r, tp and h for 7 = 0.5 as a 
function of the initial conditions. As in the text, ST and SN 
stand for Spanning Trees and Scheidegger Network respec- 
tively. SC and OCN have the same meaning as above. Values 
in parenthesis are scaling predictions. 



TQCN 



TSC 



4>OCN 



ipse 



hoc? 

ST 1.44 ±0.01 1.46 ±0.01 1.8 ±0.1 (1.79) 1.82 ±0.05 (1.85) 0.61 ± 0.03 
SN 1.44 ± 0.03 1.43 ± 0.02 1.8 ± 0.2 (1.79) 1.7 ± 0.1 (1.75) 0.61 ± 0.03 



FIG. 12. Finite size 1/L extrapolation for the value of r 
obtained with the OCN (T = 0) for the cases 7 = 0.25, 0.5. In 
both cases the initial condition is the spanning tree of Fig. [lj 
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TABLE III. Critical exponent cf> stemming from both the 
SC procedure and the OCN scheme, as a function of the value 
of 7 and of the value of the order q of ratios defined in Eqs. 
and (0) 



q 7 


= o.oo (sc ; 


)7 = o.25 (sc ; 


) 7 = 0.50 (SC ) 7 = 0.00 (OCN) 7 


1 


2.1 ±0.1 


2.2 ±0.1 


2.2 ±0.1 


2 


2.1 ±0.1 


2.2 ± 0.1 


2.2 ±0.1 


3 


2.1 ±0.1 


2.2 ±0.1 


2.3 ±0.1 



mm. 



rigues- 
Chancy, ana 



TABLE IV. Critical exponent di obtained from both the 
SC procedure and the OCN scheme, as a function of the value 
of 7 and of the value of the order q of ratios (^) and ( p"l| ) 



„ver£tty^ Press, Gjea|_]grjtain, 1997) 
[2f A. Rmaldo, I. Ro'c trigues-Iturbe, R. 



Rinaldo, Fractal River 
elf- Organization (Cambridge Uni- 



q 7 


= 0.00 (sc ; 


)7 = o.25 (sc ; 


) 7 = 0.50 (SC ) 7 = 0.00 (OCN) 7 


1 


1.3 ±0.1 


1.3 ±0.1 


1.1 ±0.1 


2 


1.3 ±0.1 


1.4 ±0.1 


1.2 ±0.1 


3 


1.3 ±0.1 


1.4 ±0.1 


1.2 ±0.1 



TABLE V. Comparison between exponents <j> as obtained 
from both the SC procedure and the OCN scheme, as a func- 
tion of the value of the initial conditions and of the value of 
the order q of ratios ( p"3| ) and (fl4|). As in the text, ST stands 
for Spanning Trees and SN for Scheidegger Network. Here 
7 = 0.5. 



q 


<j>ST (SC ) 


4>sn (sc ) 


4>st (OCN) 


0sn (OCN) 


1 


2.2 ±0.1 


1.8 ±0.1 


2.1 ± 0.1 


1.7 ±0.1 


2 


2.2 ±0.1 


1.8 ±0.1 


2.1 ± 0.1 


1.7 ±0.1 


3 


2.3 ±0.1 


1.8 ±0.1 


2.0 ±0.1 


1.7 ±0.1 



TABLE VI. Comparison between exponents di resulting 
from both the SC procedure and the OCN scheme, as a func- 
tion of the value of the initial conditions and of the value of 
the order q of ratios ( p"3| ) and (fl4|). As in the text, ST stands 
for Spanning Trees and SN for Scheidegger Network. Here 
7 = 0.5. 



q 


rfiST (SC ) 


dlSN (SC ) 


d ;sT (OCN) 


diSN (OCN) 


1 


1.1 ±0.1 


0.9 ±0.1 


1.3 ±0.1 


0.9 ±0.1 


2 


1.2 ±0.1 


1.0 ±0.1 


1.3 ±0.1 


1.0 ±0.1 


3 


1.2 ±0.1 


1.0 ±0.1 


1.3 ±0.1 


1.0 ±0.1 
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